Phase Transitions in Models for Coupled Charge-Density Waves 
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Various phase transitions in models for coupled charge- density waves are investigated by means 
of the e-expansion, mean-field theory, and Monte Carlo simulations. At zero temperature the ef- 
fective action for the system with appropriate commensurability effects is mapped onto the three- 
or four-dimensional XY model, depending on spatiotemporal fluctuations, under the corresponding 
symmetry-breaking fields. It is revealed that the three- and four-dimensional systems display a sin- 
gle transition between the clock order (with broken Zm symmetry) and disorder. The nature of the 
phase transition depends crucially on the commensurability factor M: For M > 4, in particular, the 
transition belongs to the same university class as the XY model. On the other hand, in the presence 
of misfit causing frustration in the charge-density wave, the inter-chain coupling is observed to favor 
either the commensurate state or the incommensurate state depending on the initial configuration; 
this gives rise to hysteresis around the commensurate-incommensurate transition. Boundaries sep- 
arating such phases with different symmetries are obtained in the parameter space consisting of the 
temperature, symmetry-breaking field, fluctuation strength, inter-chain coupling, and misfit. 

PACS numbers: 64.60.Cn, 71.45.Lr, 64.70.Rh 



I. INTRODUCTION 

A large number of organic and inorganic solids 
have crystalline structures in which fundamental struc- 
tural units form linear chains. 1 In these materials, 
largely different overlaps of the electronic wave functions 
in various crystallographic directions lead to strongly 
anisotropic, so-called quasi-one-dimensional (quasi-lD) 
electron bands. Among the exotic physical phenomena 
in quasi-lD materials, charge density waves (CDWs) are 
of much continual interest. Research into this topic has 
been stimulated by recent advances in experimental tech- 
niques, which now allow direct observation of CDWs and 
measurement of various static and dynamic properties of 
CDW systems^ In general the material of quasi-lD struc- 
ture is expected to exist in the form of a bundle of chains 
rather than of a single chain. In such a bundle of chains, 
inter-chain tunneling of electrons leads to coupling of the 
fluctuations on neighboring chains, which may affect the 
behavior of the system in a crucial way. For example, 
the system of two coupled incommensurate chains in the 
weak-coupling limit has been reported to exhibit a com- 
plicated commensurate-incommensurate (CI) transition, 
reminiscent of devil's staircase. 3 Obviously, the opening 
of the gap at the Fermi surface implies that each separate 
chain is an insulator at low temperatures. However, with 
the coupling between chains taken into account, expected 
are various transitions between insulating and metallic 
phases, the extensive investigation of which is still lack- 
ing in spite of the ubiquity of the CDW. 

In this paper, we investigate nature of the phase transi- 
tions in the coupled CDW system. In the absence of mis- 
fit the effective action for the (commensurate) system at 
zero temperature with suitable commensurability effects 
is mapped onto the XY model under the corresponding 



symmetry-breaking field. The effective dimension of the 
system is given by four if both spatial and temporal fluc- 
tuations are significant; otherwise it is three. The phase 
transition in the resulting three-dimensional (3D) XY 
model under the appropriate symmetry-breaking field 
is examined via the e-expansion and the Monte Carlo 
method. It is found that there emerges in the system the 
clock ordered phase (with the Zm symmetry broken) via 
a second-order transition, the nature of which depends on 
the commensurability factor M: In particular the critical 
behavior for M > 4 appears to be the same as that of 
the 3D XY model. For the four-dimensional (4D) XY 
model, on the other hand, mean-field theory is expected 
to be accurate and reveals a single transition, which is of 
the first order for M = 3 and of the second order oth- 
erwise. We then examine the effects of misfit, which not 
only change the nature of the order-disorder transition 
but also brings about a CI transition. In the presence 
of the misfit, correlations due to the inter-chain coupling 
are observed to favor, depending on the initial configura- 
tion, either the commensurate state or an incommensu- 
rate state, which gives rise to hysteresis behavior: While 
in the cooling process the incommensurate CDWs persist 
near zero temperature, in the heating process large por- 
tions of the system remain in the commensurate state at 
rather high temperatures. 

This paper is organized as follows: In Sec. [D] the ef- 
fective action at zero temperature is derived for the cou- 
pled CDW system and mapped onto the appropriate XY 
model according to whether spatial and/or temporal fluc- 
tuations are taken into consideration. Section IlIII is de- 
voted to the investigation of the phase transition in the 
3D XY model under the symmetry-breaking field, which 
introduces the Zm symmetry to the system. Here two 
independent approaches are employed: the e-expansion 
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in Sec. IIII Al and the Monte Carlo method in Sec. IIII Bl 
Section II VI presents the mean-field analysis of the 4D 
XY model, demonstrating the first- and second-order na- 
ture of the transition under the appropriate symmetry- 
breaking field. In Sec.[V]the effects of the inter-chain cou- 
pling in the presence of misfit are investigated via Monte 
Carlo simulations, which discloses properties of the CI 
transition in the system. Combining the results for the 
CI transition with those for the 3D and 4D XY mod- 
els obtained in Sec. IIIII and IIVI we construct in Sec. I VII 
schematic phase diagrams for general coupled CDW sys- 
tems in the 3D space consisting of the temperature, inter- 
chain coupling, and misfit. Finally, Sec. IVIII gives a brief 
summary. 



commensurate CDW chains reads 



II. EFFECTIVE ACTION 



We consider a system of coupled near-commensurate 
CDW chains along the z direction, each of which is 
characterized by the commensurability factor M and the 
position-dependent misfit S r with r = (x, y, z). On the 
x-y plane, the chains are assumed for simplicity to con- 
stitute a square array of lattice constant a (= 1). Dis- 
regarding amplitude fluctuations of the complex CDW 
order parameter and considering spatial and temporal 
fluctuations of the phase on\y^ we write the Hamilto- 
nian in terms of the phase </> r of the order parameter at 
position r and the momentum p r — Cdcf> r /dt: 



Hx = 



^ + ^(^T-*r) -V COS M(f> 



2 



0^ 

dz 



where for the moment the inter-chain coupling has been 
omitted. The first term and the second term correspond 
to the change of the total electron kinetic energy due to 
temporal and spatial fluctuations, respectively, whereas 
the third term represents the commensurability energy. 
The dimensionless coupling constants C, Vq and U» de- 
pend on such detailed microscopic structure of the sys- 
tem as the density of states at the Fermi level, the ef- 
fective electron mass renormalized due to the lattice vi- 
bration, the electron-phonon coupling strength, and the 
cutoff energy^ 

We now consider inter-chain tunneling of electrons be- 
tween nearest neighboring chains at (x,y) and at (x',y') 
on the x-y plane, i.e., at the same position z; this gives 
rise to the interaction of the form U±_ cos(0 r — </) r >), where 
U± is the dimensionless inter-chain coupling constant 
and higher-order harmonics have been disregarded. With 
this included, the Hamiltonian TL for the coupled near- 



n = 




■Br + iri-Er— 5 *} - cos M<j), 



u± cos( 



(2) 



(xy,x'y') 



where the position <j) r and the conjugate momentum 
p r are considered to observe the commutation relation 
[4>r,Pr>] = «<5r,r', suggesting the position representation 
p r = — id/d(j> r . In the second summation (xy,x'y') stands 
for the nearest neighbor pairs on the x-y plane at fixed 
z. Throughout this work, we set h = 1, c = 1, and the 
Boltzmann constant kg = 1 . 

Note that in the presence of the commensurability en- 
ergy, the misfit cannot be simply gauged away and in- 
troduces frustration to the system. For the time be- 
ing we consider the case of strictly commensurate CDW 
systems without any misfit. To investigate the quan- 
tum phase transition in this case, driven by quantum 
fluctuations at zero temperature, we follow the standard 
procedure^ to map a c?-dimensional quantum system to 
a (d+l)-dimcnsional classical system and obtain the cor- 
responding effective action. In Sec. Ill Al the system with 
negligible spatial fluctuations along the chain direction is 
considered and the corresponding effective action at zero 
temperature is mapped onto the 3D XY model, where 
commensurability effects are described by the appropri- 
ate symmetry-breaking field. Section III Bl discusses the 
system in the presence of spatial fluctuations. We first 
pay attention to the classical limit, where temporal fluc- 
tuations may be neglected. Here the system is intrinsi- 
cally 3D, and the effective action is again mapped onto 
the 3D XY model. Next, with both spatial and temporal 
fluctuations considered, the appropriate effective action 
at zero temperature is identified with the 4D XY model. 



A. Homogeneous Case 

Although it is in general expected that strong fluctua- 
tions are present, we for the moment assume that spatial 
fluctuations along each chain are negligible, which gives a 
two-dimensional (2D) system without z dependence. In 
the absence of the symmetry-breaking field due to com- 
mensurability, Eq. J5J reduces to the 2D XY model with 
kinetic energy 



n 



^ 2C 



E 

(r,r'> 



U± cos( 



'). 



(3) 



where the summation in the second term is to be per- 
formed over all nearest-neighboring pairs in the 2D space 
with r = (x,y). The Hamiltonian in Eq. © has been 
studied in the context of quantum arrays of Josephson 
junctions. 6 In particular, at zero temperature the 2D 
quantum system in Eq. (0 is well known to map onto 
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a 3D classical system via the standard lore^ Introducing 
the imaginary time r axis and dividing the interval be- 
tween r = and r = T _1 into N slices of equal width 
At = 1/NT, in the zero temperature limit (T — » 0) we 
arrive at the partition function of the anisotropic 3D XY 
model 



T>4> exp 



5> r 

<r,r'} 



(4) 



where r = (r, x, y) represents the position in the 3D 
space, consisting of the (imaginary) time r and the 2D 
space (x,y). The anisotropic coupling is defined on each 
bond: 



r,r' 



C/Ar for r' = r ± f At, 
U±At for r' = r ± x or r ± y. 



Strictly speaking, we should keep Ar infinitesimal. With- 
out affecting the universality, however, we can rescale the 
space and time and obtain the partition function of an 
isotropic 3D XY model 



Z = <j)Vcj)e- H 



(5) 



with the desired effective action 



— H = K 2_, cos (0r — <j>r>)i 
<r,r') 

where K~ x = (Ct/i) -1 / 2 measures the amount of quan- 
tum fluctuations. 

We then accommodate the commensurability effects at 
zero temperature, taking advantage of the fact that the 
symmetry-breaking field term is diagonal in the position 
basis. Via the same procedure and rescaling as in the ab- 
sence of the symmetry-breaking field, the effective action 
for the 3D XY model obtains 



-H 



K cos( 

<r,r'> 



') + h^2cosM<p r , (6) 



where the symmetry breaking field is given by h 
V At = V 0V ^/UI- 



B. General Case 

We now turn to the system in which spatial fluctua- 
tions in the z direction are present (but still without the 
misfit), and first consider the case of only spatial fluctu- 
ations, with temporal fluctuations negligible. This corre- 
sponds to the classical limit in the sense that the momen- 
tum and the position decouple and the momentum part, 
which can be integrated out, does not affect the relevant 
physics. In this intrinsically 3D case, it is revealing to 
resort to the discrete formulation and replace the inte- 
gration J dz by the summation ^ 2 Az with sufficiently 



small Az. Regarding spatial fluctuation-dependent en- 
ergy as the continuum form of the cosine interaction in 
the discrete representation, we obtain from Eq. J5J the 
effective Hamiltonian for the 3D XY model: 



TL = Kr,r> cos(0 r — <f) r i 

<r,r') 



cos Mcj> T , (7) 



where r = (x, y, z) represents lattice sites in the 3D space, 
and the coupling strength and the symmetry-breaking 
field are given by 



{U±Az for r' = r ± x or r ± y 
U\\/Az for r' = r ± zAz 
h = V Az. 



Interestingly enough, Eqs. © and Q) show that both the 
(zero-temperature) quantum phase transition in the ab- 
sence of spatial fluctuations and the (finite-temperature) 
classical phase transition in the absence of temporal fluc- 
tuations are described by the same 3D XY model under 
symmetry breaking fields. Quantum fluctuations in the 
former play the role of thermal fluctuations in the latter, 
tending to restore symmetry. 

We finally consider the general case, where both spatial 
and temporal fluctuations are significant. Adding the 
kinetic energy term to the 3D effective action in Eq. (JJJ) 
and following the same procedure as that in Sec. Ill Al at 
zero temperature, we obtain the 4D XY model with the 
effective action 

-H = K it j' cos(<f>i - ) + h cos M(j)j, (8) 

(i,i>) i 

where i = (r, x, y, z) denotes 4D space-time lattice sites 
and the anisotropic coupling and the symmetry-breaking 
field are given by 



Az 
C 



for i' = i ± iAz 



for i' = i ± tAt 



AzAt 

AzAtU± for other nearest neighbors 
VqAzAt 



In this manner, the coupled near-commensurate CDW 
systems (without misfit) can be described by the appro- 
priate XY models under symmetry-breaking fields. De- 
pending on whether spatial and/or temporal fluctuations 
are present, the effective dimension of the system is de- 
termined to be three or four: It is four if both fluctua- 
tions are significant at zero temperature and three oth- 
erwise. The commensurability effects are described by 
the symmetry-breaking field, which introduces Zm sym- 
metry to the system. Such a symmetry-breaking field 
affects the ground-state symmetry and is expected to be 
relevant in the sense of the renormalization group (RG) 
theory. In the limit M — > oo, however, the Zm symmetry 
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is hardly distinguishable from the underlying U(l) sym- 
metry in the XY model. Therefore we expect the phase 
transition to be crucially dependent upon the commen- 
surability factor M, and devote the next two sections to 
the investigation of the phase transitions in the 3D and 
4D XY models under symmetry-breaking fields. 



with the weight factor (u > 0) 



exp[-W(S)] = exp 



uS 4 



III. 3D XY MODEL UNDER 
SYMMETRY-BREAKING FIELD 

In this section we investigate the phase transition in 
the 3D XY model whose effective action is given by 
Eq. JBJ. In the absence of the symmetry-breaking field 
(ft = 0), the 3D XY model has been studied both 
analytically^^ and numerically, revealing that vortex 
loops do exist and proliferate at the phase transition. 
Accordingly, the topological scaling idea has been ex- 
tended to the 3D transition with conventional long-range 
order, with the scaling procedure of Ref. ll] generalized 
appropriately for 3D directed loops&S One may then be 
tempted to extend the study of the 2D model in Ref. [l2j 
to incorporate the symmetry-breaking term in the 3D XY 
model, and combine recursion relations for the vortex fu- 
gacity ?/o & n d the field perturbation according to the 
duality relation between them. However, in contrast to 
the 2D case, the geometric scaling of the coupling con- 
stant in three dimensions results in nonzero values of yo 
and yh at the fixed points 

On the other hand, the e-expansioni^ is well known to 
provide a mathematical formalism for calculating critical 
exponents of the O(n) model near four spatial dimen- 
sions, allowing classification of the universality class of 
the system. Here we extend the e-expansion approach to 
incorporate the effects of the symmetry-breaking field. 
The e-expansion turns out to be useful only for small 
commensurability factor M; we thus supplement its lim- 
itations with the Monte Carlo numerical method. In 
Sec. lIII Al the symmetry-breaking perturbation is treated 
within the e-expansion and Sec. 1111 Bl is devoted to the 
Monte Carlo simulations of the 3D XY model under 
symmetry-breaking fields. 



A. e- Expansion 



which has been expanded up to 0(S 4 ). At the physical 
XY fixed point the sixth-order term is less relevant in 
dimension d = 4 — e than d = 4>i4 With the identification 
= cos 4> T and Sj! — sin <f> T , the original action in Eq. © 
is written in terms of the continuous spin S r : 



h = x^-s^x 

r,r' r 

M(M-3) 



M 

{2S*) M ~(2S*) M - 2 



2! 



(2fi 



x\M-4 



(10) 



Note that for M < 4 the order of the action does not 
exceed 0(S 4 ), leading to the Ginzburg-Landau- Wilson 
(GLW) effective Hamiltonian H = H + W to 0(S 4 ). 

We first examine the case M = 2. In the momen- 
tum space representation, keeping only relevant terms to 
0(q 2 ) and scaling the spin variable according to S q = 
(Aa d+2 )~ 1//2 er q , we express the GLW Hamiltonian as 



-H = 



+ Vi f f I <<<"^ q - q - q " (n) 

J q J a' J a" 



+ 2V 2 



+ V 3 



qJq'Jq 



^q'^-q-q'-q'" 



where the coefficients are given by 



We consider the e-expansion for the 3D XY model, ex- 
tending the original formulation for the O(n) modeli^ to 
incorporate the symmetry-breaking field. To begin with, 
we employ the two-component continuous local variable 
or spin S r = (S*,S}!) at each 3D lattice site r, the total 
number of which is denoted by N. The constraint that 
each spin has the unit magnitude is relaxed by the addi- 
tional weight factor introduced to the partition function: 



n 



dS r 



-[H+W(S r )] 



(9) 



1 



Ka 2 
1 

AV 



(b-dK- 2ft) 



(b - dK) 



V! = V 2 = v 3 



-u<0 



and J q = / dq/(2-K) d = (1/Aa rf )£ q with the lattice 
constant a restored for clarity. Here it is shown that the 
M = 2 symmetry-breaking field gives rise to anisotropy 
in the quadratic term, making r x less than r y . Following 
the standard procedure, we obtain the recursion relations 
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to the leading order in e = 4 — d: 



dr x 
dl 
dry 
dl 
dVi 
~df 
dV 2 
dl 



2rv 



12C 
1 + rv 



Vi 



4C 
1 + r v 



Vn 



4C Tr 12C Tr 
2r y - - V 2 - — V 3 



eV 1 
eV 2 



l + r x 
36C 

16C 



l + r v 



V, 2 



4C 



-V 2 

2 K 2 



(1 + ^(1 + r„) 
12C 



(l + r x )(H-r v ) 



W 2 + 



12C 



(l + r x )(H-r v ) 



(12) 



W 3 

9; 



-ST = ^3 



36C 



4C 



{l + r x f 3 T (l + r,)2^ 2 



with the spatial scale factor I and C = 
21-^/2/(^/2-1)!. 

In the absence of the symmetry-breaking field (h = 0), 
the parameter space reduces to the 2D space (r, V") since 
r 2 , = r y = r and Vi = V2 = V 3 = V. In this case the 
nontrivial fixed point of the recursion relation in Eq. I|12|) 
is simply the XY fixed point, given by 



■5' 



V* 



40C 



(13) 



to 0(e). As the symmetry-breaking field is turned on, 
however, we have r x < r y at the initial locus and need to 
examine the RG flow in the full five-dimensional (5D) pa- 
rameter space (r x ,r y , V\, V2, V 3 ). To investigate the sta- 
bility of the XY fixed point in the 5D parameter space, 
we linearize Eq. 11 2[) about this point and obtain a sta- 
bility matrix which yields five eigenvalues together with 
corresponding scaling exponents: 



yi 



e 

5' 

V4 " 



2c 
J 1 



Ac 



(14) 



■5' 2/5 = ~T 



Since the exponents 2/3, j/4, and 2/5, which correspond 
to the three eigenvectors spanning the 3D subspace 
(Vi, V2, V3), are all strictly negative, it is concluded that 
all Vi's (for i = 1,2,3) are irrelevant at the XY fixed 
point. In contrast, both y\ and j/2 are positive, indicating 
that the XY fixed point is unstable in the 2D subspace 
(r x , r y ). The initial locus, which is kept off the XY fixed 
point by the symmetry-breaking field, should flow far off 
the XY fixed point. Instead it is expected to flow toward 
the Ising fixed point located at 



e 

T 

6' 



v 



00, v? = -—, v 2 * = v 3 * = o. 

(15) 

At this point it is more appropriate to take r y 1 as the 
scaling field, giving the scaling exponents 



c 

3' 



-2. 



(16) 



Thus the initial locus with r x < r y indeed flows toward 
the Ising fixed point along the stable r y direction. In 
this manner, the symmetry-breaking field for M — 2 in- 
troduces anisotropy in the quadratic term, making the 
XY fixed point unstable and generating an Ising or two- 
state clock fixed point. Revealed accordingly is a sin- 
gle second-order transition between two-state clock order 
and disorder. 

We next turn to the case M = 3, where Eq. lflTj|) shows 
that terms of linear and third order in the x-component 
of the spin come into play; power counting suggests that 
these fields are relevant near d = 4. Owing to the 
anisotropy associated with the absence of S v terms, in 
particular, the cubic term (S x ) 3 here may not be re- 
moved by mere shift and is not redundant, in contrast 
to the case of the Ising model. The symmetry-breaking 
field is thus expected to drive the transition between the 
disordered phase and the three-state clock ordered phase. 

Finally the symmetry-breaking field for M — 4 intro- 
duces anisotropy in the quartic term: 

V x = -u + 8a d - 4 K- 2 h 
V 2 = V 3 = -u. 



(17) 



This leads the action in Eq. (|11|) to be unstable for suffi- 
ciently large values of the field h, making it necessary to 
consider higher-order terms in the weight function W(S) 
for stability. Unfortunately, such calculation of higher 
cumulants is quite a formidable job, and it is very dif- 
ficult to extend the e-expansion to the case of a high 
commensurability factor. 



B. Monte Carlo Simulations 

This section presents the Monte Carlo study of the 3D 
XY model in symmetry-breaking fields. To estimate the 
critical temperatures and determine critical exponents, 
we have performed Monte Carlo simulations at "temper- 
atures" (i.e., quantum or temporal fluctuations) ranging 
from K^ 1 = 0.5 to K^ 1 — 5 on lattices of linear size 
L = 4 up to L = 32 for several commensurability factors 
and field strengths. Measured in simulations are the or- 
der parameter m and the susceptibility \ defined to be 



X 




N ^ 



,i<t>t 



(18) 



where TV is the number of sites. We have employed the 
multiple histogram method^ to interpolate the quanti- 
ties calculated sparsely in a given range of the temper- 
ature to any temperature inside the range. It not only 
saves a great deal of computing time but also provides the 
values of quantities at arbitrary temperatures for finite- 
size scaling, resulting in critical exponents of better ac- 
curacy. In fact we have obtained the best collapse of the 
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h 



FIG. 1: Data collapse of the order parameter m for three 
different values of the system size L, with the best-estimated 
critical temperature and exponents given in the text. The 
inset shows the dependence of the error on the value of v, 
with T c and fixed at their best estimated values. The data 
have been taken for M — 2 and at h — 2.7; typical error bars 
are not larger than the symbol sizes. 



scaling function such as 

m(j}l v t) = L p ' v m{t) (19) 

with the reduced temperature t = T/T c — 1 = K c /K — 1, 
by minimizing the measure of error— 



1 



dx 



«m 2 (x))) - ((m(x))Y 



( 2 °) 

over the critical temperature and exponents. Here ((■ • • )) 
stands for the average over the lattice size and similar 
relations for the susceptibility have also been considered. 
We have taken the range [x m i n ,x max ] = [—Ax, As] and 
attempted a series of collapses as Ax is diminished. The 
result has then been extrapolated to the limit Ax — > 0. 

Figure ^ presents typical collapse of the scaling func- 
tion with the best-estimated critical temperature T c and 
exponents and v\ nice collapse behavior can be ob- 
served near the critical temperature (t — 0). The inset 
in Fig. ^ discloses how the corresponding error cr^ de- 
pends on the value of the exponent v. We have thus es- 
timated the error in the obtained critical exponent from 
the standard deviation of the minimizing values over the 
sample ensemble. In this way, the phase transition for 
71/ = 2 is found to be of the second order with the crit- 
ical exponents v = 0.63 ± 0.01 and = 0.34 ± 0.04 for 
h = 1.2 and v = 0.627 ± 0.004 and = 0.32 ± 0.01 for 
h = 2.7. These results coincide perfectly with the known 
critical exponents for the 3D Ising model: v — 0.630 
and — 0.324, thus demonstrating the validity of the 
e-expansion analysis in Sec. IIII Al Here it is of interest 
that the same universality class was also reported in the 
quantum Monte Carlo study of a tight-binding model of 
spinless fermion chains coupled by intra- and inter-chain 
Coulomb interactions— The tight-binding model anal- 



FIG. 2: Critical temperature K~ x in the 3D XY model under 
symmetry-breaking fields versus the symmetry-breaking field 
strength h, for the commensurability factor M — 2,3, 4, and 
5. Typical error bars are smaller than the symbol size and 
lines are merely guides to eyes. 



ysis, though being restricted to the M = 2 case, could 
incorporate amplitude fluctuations of the CDW order pa- 
rameter directly. Indeed the accordance between the two 
analyses manifests that fluctuations in the phase (rather 
than in the amplitude) mostly determine the nature of 
the transition. 

For M = 3, on the other hand, we have found no 
evidence for the first-order transition, and obtained the 
estimation v = 0.596 ± 0.007 and = 0.3 ± 0.01; this is 
to be compared with the long-standing controversy as to 
whether the transition in the 3D three-state Potts model 
is of the first order or continuous— It is also of interest 
that the critical behavior for M > 4 is similar to that of 
the 3D XY model without any symmetry-breaking field. 
Figure El shows the critical temperature (i.e., the critical 
fluctuation strength) K c versus the field strength h for 
2 < M < 5. It is observed that for M < 3 the critical 
temperature increases with the field strength while for 
M > 4 the commensurability energy appears to have no 
effect on the critical temperature, in the range probed. 
Such critical behavior for M > 4, similar to that of the 
ordinary 3D XY model, reflects that, as noted in Sec. till 
the Zm symmetry for large M is indiscernible from the 
U(l) symmetry underlying in the XY model. 



IV. 4D XY MODEL UNDER 
SYMMETRY-BREAKING FIELD 

We now study the phase transition in the 4D XY 
model, described by the effective action in Eq. |JSJ, 
through the use of the mean-field approximation. Here 
the mean-field approximation, developed for supercon- 
ducting arrays in applied magnetic fields^ is expected 
to be accurate since the upper critical dimension of the 
XY model is given by d u — 4. With the space and time 
rescaled appropriately, the (mean-field) self-consistent 
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equation reads 



271 



(21) 



where Zmf is the partition function 

,-2tt 







corresponding to the mean-field action 

—Hmf = 4JsT((cos0) cos</> + (smcj))) smcj) + hcosMtfi 

(22) 

with the rescaled coupling constant K = \fCU\_ and field 

h = v ^cjul. 

In the absence of the symmetry-breaking field, the self- 
consistent equation leads to the equation of state for the 

order parameter m = \/ (costfi) 2 + (smcj)) 2 : 



Iu{AKmY 



(23) 



where /„ is the nth modified Bessel function. The system 
is in the disordered phase for K < K c (h=0) = 1/2 char- 
acterized by m = 0; beyond K c (h=0) the ordered phase 
with m ^ is favored. Similarly, the M = 1 case can 
be analyzed easily, where arbitrarily small but positive h 
results in (cos0) > and (sin</>) = 0. The equation of 
state is thus 



h(AK (cos^) +h) 
I {AK (cos0) +h)' 



(24) 



For larger values of M, Eq. (|21|l may be solved nu- 
merically with the parameters K and h varied. It is 
found that there exists a field-dependent transition cou- 
pling strength K c (h) below which Eq. H21|) bears only 
the null solution (cost/)) = (sin</>) =0. As K is increased 
beyond K c {h), nonzero stable solutions emerge. Due to 
the Zm symmetry, i.e., the invariance of the action un- 
der the shift in the angle by 2-kti/M for any integer n, 
the solution of Eq. I|21[l can be expressed in terms of the 
M-state clock order parameter m(K, h): 



i<f>\ — 



<e*> 



m(K, h)e 2m/M 



(25) 



with integer n = Q, 1, . . . , M — 1. The critical coupling 
strength K c (h) is defined to be the largest value of K 
satisfying m(K, h) = 0. 

Figure [21 shows how the order parameter m depends 
on the parameters K and h for different values of M. It 
is observed that m in general increases rapidly from zero 
as K exceeds K c and then saturates to unity, with the 
increase more rapid for larger h. While for Al = 2 and 
M > 4 the transition is found to be of the second order, 
numerical results for M = 3 indicate a first-order transi- 
tion, which is expected from the appearance of the third- 
order term in the e-expansion. For M < 3, the critical 
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FIG. 3: Order parameter m as a function of the coupling 
K and the field strength h for (a) M = 2, (b) M = 3, 
and (c) M — 4. Main graphs show m versus K for 
h = 0(+),0.5(x), 1(D), 5(0), and 10(A) from below. In- 
sets display m(K, h) versus h for various values of K: From 
below, (a) K = 0.275,0.3,0.375,0.5, and 0.75; (b) K = 
0.47, 0.475, 0.4875, 0.5, and 0.75. The overall behavior of the 
order parameter for M > 4 is the same as that for M = 4. 
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FIG. 4: Phase diagram for the 4D XY model under the 
symmetry-breaking field, displaying the boundaries between 
the M-state clock order phase (M) and the disordered phase 
(D) for different values of M. Lines for M = 2 and 3 are 
merely guides to eyes. 



coupling strength K c (h), starting from K c (h=0) = 1/2, 
decreases with h almost exponentially to the asymptotic 
values, 0.250 and 0.462 for M = 2 and for M = 3, re- 
spectively. For M > 4, on the other hand, the transition 
point does not depend on the strength of the symmetry- 
breaking field, giving K c = 1/2 regardless of h. Again 
manifested is the crucial role of commensurability in the 
phase transition. 

The resulting phase diagram for the 4D XY model 
under the symmetry-breaking field is displayed in Fig. 01 
where the boundary K~ versus h is plotted for different 
values of M. When K~ x < K~ l , the Zm symmetry as 
well as the U(l) symmetry is broken and identified is the 
M-state clock ordered phase, with one of the minima of 
—hcosMcj) favored. When K^ 1 > K~ x , unlike the U(l) 
symmetry broken explicitly for h ^ 0, the Zm symmetry 
remains unbroken, leading to the disordered phase. 



V. COMMENSURATE-INCOMMENSURATE 
TRANSITION 

Up to the present, we have assumed the absence of 
misfit, so that only the transition between the disor- 
dered phase and the M-state clock ordered phase, where 
commensurate CDWs are developed, has been consid- 
ered in several idealized cases described in Sec. [H] How- 
ever, the misfit, being a key ingredient to bring about 
the commensurate-incommensurate transition in a near- 
commensurate CDW chain, must be taken into consider- 
ation in understanding various transitions in the coupled 
CDW system. 

It is well known from the study of the ID Frcnkcl- 
Kontorova modeli^ at zero temperature that a single 
near-commensurate CDW changes from the commensu- 
rate state to the incommensurate state when the misfit 
exceeds a critical value depending on the commensurabil- 



ity energy. At finite temperatures the ID CDW system 
is always in the incommensurate stated On the other 
hand, in the coupled CDW system, interactions between 
the CDW chains may alter the nature of the transition. 
First of all, the effective dimensions of the system grow 
to three, giving rise to the persistence of long-range or- 
der, as observed in the previous sections. Moreover, the 
CI transition itself can also be affected by the inter-chain 
coupling in that the interactions favor either the commen- 
surate state or the incommensurate state, as explained 
below. To investigate the phase transition of the coupled 
CDW system, we consider the Hamiltonian in Eq. (J2J 
without temporal fluctuations, which makes the problem 
simpler and allows to focus on the static properties only. 

At zero temperature the problem is rather simple to 
solve. The inter-chain coupling term in Eq. © reaches 
the minimum when all the phases <f> T in the xy-plane be- 
come equal to each other. All the CDW chains, therefore, 
follow the same phase configurations determined by the 
ID Frenkel-Kontorova model, thus undergo the CI tran- 
sition at the critical misfit given by^ 



(26) 



regardless of the inter-chain coupling strength U±. 

The system at finite temperatures is examined by 
means of the Monte Carlo method. For this we discretize 
the z-axis in Eq. as in Sec. Ill Bl and obtain the lattice 
Hamiltonian 

H = ^^r+z-^r-^-^VuCOSM^. 




H Y U±COs(4> T -(f)r>), 
z (xy,x'y') 



(27) 



where Az has been absorbed into the coupling constants. 
Since the dimension along the CDW chain is usually 
longer than the lateral dimensions, the lattice size L z 
along the z-axis is kept to be two times larger than the 
other linear sizes in our simulations, and accordingly the 
inter-chain coupling U±_ is restricted to be smaller than 
Un to avoid inessential finite-size effects. We further set 
the parameters to U\\ = 1 and Vq — 0.2 and sweep the 
inter-chain coupling strength from 0.01 to 0.5 and the 
misfit 5, which is assumed to be uniform throughout 
the system, from zero to 0.8 beyond the critical misfit 
5 C (= 0.569). 

The order-disorder transition is described conveniently 
by the (incommensurate) order parameter defined to be 



N t— 1 



i{4> r —6z) 



(28) 



and its susceptibility x- Note that in the thermodynamic 
limit this order parameter vanishes not only in the dis- 
ordered phase but also in the commensurate phase; in 
the system of finite size it may remain non-zero even 
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FIG. 5: Behavior of the order-disorder transition tempera- 
ture Td with the misfit 8 in the M = 2 coupled CDW sys- 
tem for U\\ = 1 and Vb = 0.2. Each symbol corresponds 
to a different value of the inter-chain coupling strength: 
U± = 0.01(D), 0.02(B), 0.05(o), 0.1(«), 0.2(A), 0.3(A), 0.4(v), 
and 0.5(T). Typical error bars are smaller than the symbol 
sizes and lines are merely guides to eyes. 



in the perfectly commensurate state. We also compute 
the soliton density p s r^ i.e., the soliton number per site 
per chain, which characterizes the CI transition. For the 
detection of any hysteresis, these two physical quanti- 
ties have been measured in two different ways: in the 
cooling-down process with randomized initial phase con- 
figurations and in the heating-up process starting from 
the zero-temperature ground state. In each process the 
system has been equilibrated sufficiently while the tem- 
perature is varied gradually. 

Shown in Figgis the dependence of the order-disorder 
transition temperature Td on the misfit 5 in the M = 2 
CDW system at various inter-chain coupling strengths. It 
is observed that the misfit tends to reduce the transition 
temperature: The transition temperature first decreases 
as the misfit is increased from zero, then saturates when 
the misfit reaches a value comparable to S c . Here stronger 
inter-chain coupling, which helps to increase the transi- 
tion temperature, in general weakens the effects of the 
misfit. Besides the transition temperatures, the critical 
exponents also change with the misfit. For instance, the 
critical exponents in the system with U± = 0.1 are found 
to be v = 0.79 ± 0.02 and 7 = 1.33 ± 0.01 at 5 = 0.2; 
v = 0.77 ± 0.03 and 7 = 1.64 ± 0.05 at 8 = 0.4. This 
suggests that the introduction of the misfit changes the 
nature of the ordered phase, making it different from the 
3D XY ordered phase. It should be identified as the 
incommensurate CDW phase, as demonstrated below. 

We now draw our attention to the CI transition oc- 
curring at lower temperatures. In the absence of misfit 
(S = 0) the soliton density p s is observed to vanish on 
average at all temperatures below Td , indicating that the 
CDWs formed are commensurate even in the presence of 
thermal fluctuations. Such a commensurate phase is de- 
stroyed by the introduction of misfit, even for 6 < 5 C , with 
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FIG. 6: (a) The order parameter m and (b) the soliton density 
p s versus the temperature T in the M = 2 system of misfit 
8 = 0.2 and size L z — 32. Each symbol corresponds to a 
different value of the inter-chain coupling strength as listed in 
Fig. 1^] The inset is an enlarged view of the soliton density in 
the range 0.1 < T < 0.5. 



the help of thermal fluctuations. Figure[S]exhibits the be- 
haviors of the order parameter and of the soliton density 
as the temperature is varied in the system with misfit 
S = 0.2. It is shown that the order parameter reaches its 
maximum at temperature near Tjj and decreases to zero 
as the temperature is reduced while the soliton density 
begins to decrease from its maximum value at T ps Td 
and vanishes almost to zero at very low temperatures. 
Nonzero values of the order parameter and of the soliton 
density together in the regime T < To correspond to the 
formation of incommensurate CDWs. Thus the order- 
disorder transition at T = To is identified as the tran- 
sition between the incommensurate CDW state and the 
disordered state. In addition, the vanishing soliton den- 
sity [see the inset in Fig. |Sj is a sign of the commensurate 
CDW phase; this indicates the presence of a CI transition 
driven by thermal fluctuations. The CI transition tem- 
perature Tic is defined to be the temperature at which 
the soliton density begins to be non-zero (in the thermo- 
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dynamic limit). Although our computing ability disallow 
us to determine the precise value of Tjc, Fig.|SJb) man- 
ifests that Tjc increases with the inter-chain coupling 
strength. In this sense the inter-chain interactions favor 
the commensurate state. In contrast, it is also observed 
that the inter-chain interactions prefer the incommensu- 
rate state at the misfit values somewhat larger (but still 
smaller than 5 C ). 

Figure |7| compares the order parameter and the soliton 
density at 5 = 0.4, measured in two different processes: 
the cooling process and the heating one. Unlike the case 
S = 0.2, a hysteresis is evident at 5 = 0.4 (and also in the 
undisplayed case of 6 = 0.3), as shown in Fig. [3a) and 
(b) for the order parameter and in Fig. (Jfc) and (d) for 
the soliton density. In the heating process staring from 
the commensurate ground state at zero temperature, the 
commensurate CDWs survive some thermal fluctuations 
and experience the CI transition at a finite transition 
temperature. In contrast, as the temperature is lowered 
in the cooling process, the order parameter keeps increas- 
ing and the soliton density saturates to a finite value, un- 
less the inter-chain coupling is sufficiently weak (i.e., for 
U±_ > 0.02). This implies that the system still consists of 
the incommensurate CDWs at temperatures where com- 
mensurate CDWs are supposed to be more stable. For 
weak coupling (U± < 0.02), on the other hand, no dis- 
crepancy between the two processes is observed. 

It is of interest to compare the hysteresis observed here 
with the one reported in the specific heat around the CI 
transition in a number of incommensurate systems^ The 
latter is attributed to the effects of pinning of the incom- 
mensurate modulation due to defects or impurities; our 
observation, in contrast, shows that the hysteresis can ap- 
pear even in the absence of defects, for which inter-chain 
interactions are responsible. Namely, the inter-chain cou- 
pling operates in different ways depending on the pro- 
cess: In the cooling process, correlations between CDW 
chains due to the inter-chain interactions hinder each 
CDW chain from getting into the commensurate state. 
On the contrary, in the heating process, the inter-chain 
correlations hold each chain close to the zero-temperature 
ground state until thermal fluctuations become compara- 
ble to the inter-chain interaction energy. This argument 
is supported by the observation that the order parame- 
ter and the soliton density increase abruptly in a narrow 
region of the temperature, as shown in Fig.Qb) and (d), 
and such discrepancy becomes manifest as the inter-chain 
coupling becomes stronger. In addition, in the same nar- 
row region the soliton number and the order parameter 
fluctuate strongly while the energy fluctuations are rela- 
tively weak. This suggests a kind of configurational fluc- 
tuations over many metastable states. 

For large values of the misfit (S > S c ), neither com- 
mensurate CDW nor hysteresis is observed and the soli- 
ton density does not vanish even at zero temperature in 
any process. Only the transition separating the incom- 
mensurate CDW state from the disordered state is thus 
identified. 



We have also performed similar simulations for the 
commensurability factor M = 3, only to obtain results 
qualitatively the same as those for M = 2 presented 
above. Quantitative differences observed include that for 
M = 3 the misfit affects the transition temperature more 
weakly and the hysteresis takes place even at smaller mis- 
fit such as 5 = 0.2. 



VI. PHASE DIAGRAM 

Combining the results for the 3D and 4D XY models 
under symmetry-breaking fields, studied in Sees. ITTll and 
IIVI and those for the near-commensurate CDW model, 
studied in Sec. [V] we are ready to describe the phase 
transition in the coupled CDW system. First, we fo- 
cus on the system in the absence of any misfit (6 = 0). 
For this purpose, it is adequate to consider the 3D space 
consisting of the three parameters: the temperature T, 
the strength of the symmetry-breaking field h, and the 
amount of quantum or temporal fluctuations A -1 . On 
each of the three planes (T = 0, A -1 = 0, and h = 0) in 
the 3D space, appropriate boundaries separating various 
phases can be plotted through the use of the known re- 
sults, as shown in Fig. [S] The phase at a given point in 
the 3D space can then be speculated by means of finite- 
size scaling in the region near T — or the semi-classical 
methods in the A _1 w region. 

We first consider the homogeneous case, where spatial 
fluctuations along each chain are negligible, and show 
the phase diagram in Fig. [S] for (a) M = 4 and (b) 
M > 5. At zero temperature, the system is mapped onto 
the 3D XY model as discussed in Sec. |n] and the corre- 
sponding phase diagram, obtained in Sec. lIIII is sketched 
on the h-K~ x plane. On the other hand, in the ab- 
sence of temporal fluctuations (AT -1 = 0), the system 
maps onto the classical 2D XY model under symmetry- 
breaking fields, for which the phase diagram in Ref. [l2j 
is drawn on the h-T plane. Finally, on the T-K^ 1 plane 
with h — 0, the system is described by the Hamilto- 
nian in Eq. J2J and expected to display the Berezinskii- 
Kosterlitz-Thouless transition renormalized by quantum 
fluctuationsiSiS^ Note that the cases M = 2 and 3 are not 
shown here. In this case presumably Ising/Potts critical 
lines exist on the h-T plane; however, it is not known how 
these lines connect up to the phase boundary for h = 0. 

In the general case with spatial fluctuations present, 
we obtain the phase diagram shown in Fig. Etc) and (d) 
for M < 3 and M > 4, respectively. Owing to the ad- 
ditional dimension along the chain direction, the system 
maps onto the 3D XY model in the classical limit, i.e., 
on the h-T plane, whereas 4D XY model is obtained at 
zero temperature. Accordingly, on the h-K^ 1 plane, the 
phase diagram of the 4D XY model obtained in Sec. IIVI 
is drawn. It is observed that the 3D XY ordered phase 
does not survive the field and is replaced by the M-state 
clock ordered phase for h ^ 0. For M > 4, however, 
the commensurability energy affects neither the transi- 
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FIG. 7: The order parameter m [(a) and (b)] and the soliton density p 3 [(c) and (d)] versus the temperature T in the M — 2 
system of misfit S = 0.4 and size L z = 32. The left figures [(a) and (c)] and the right ones [(b) and (d)] show the data obtained 
in the cooling process and in the heating process, respectively. Each symbol corresponds to a different value of the inter-chain 
coupling strength as listed in Fig. 



tion temperature nor the nature of the transition. 

We next take into consideration the effects of misfit 
and show in Fig.[!|]the schematic phase diagram in the 3D 
space consisting of T, UJ_ , and S, with U\\ and Vo fixed. 
At 6 = 0, the system belongs to the same university class 
as the 3D XY model under the symmetry-breaking field, 
for which the phase diagram is drawn on the T — U^ 1 
plane. On the other hand, at zero temperature, each 
CDW chain undergoes the CI transition with the critical 
misfit S c which is independent of the inter-chain coupling 
strength U±. The order-disorder transition temperature 
To and the CI transition temperature Tjq obtained in 
Sec. [V] produce surfaces depicted by (thick) solid and 
dashed lines, respectively, in the (T, ,<J) space. The 
phase diagram shows that for 6 < 8 C the system under- 
goes double transitions as the temperature is lowered: 
first, the order-disorder transition into the incommen- 
surate CDW phase and then the CI transition into the 
commensurate CDW phase. As revealed in Sec.0 a hys- 
teresis takes place around the CI transition point, which 



manifest itself more clearly as the inter-chain coupling 
becomes stronger. 



VII. CONCLUSION 

To study phase transitions in coupled CDW systems, 
we have mapped the systems at zero temperature onto 
three- or four dimensional XY models, depending on the 
spatiotemporal fluctuations, under symmetry-breaking 
fields which arise from the commensurability energy. 
Such techniques as e-expansion, mean-field theory, and 
the Monte Carlo method have then been applied to the 
obtained XY models. Revealed is a single second-order 
transition between the M-state clock order and disorder 
in both the three-dimensional and four-dimensional sys- 
tems, except for the case M = 3 in the four-dimensional 
system where the transition is of the first order. In par- 
ticular, the commensurability with M > 4 has been ob- 
served not to change the properties of the transition in- 
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FIG. 8: Schematic phase diagram for the coupled commensurate CDW system in the {T,h,K~ l ) space, depending on spa- 
tiotemporal fluctuations and commensurability: the homogeneous case (a) for M = 4 and (b) for M > 5; the general case (c) 
for M — 2 or 3 and (d) for M > 4. Plotted on each plane are boundaries between various phases, including the disordered 
phase (D), the M-state clock ordered phase (M), the algebraically ordered phase present in the 2D XY model (2D XY), and 
the 3D XY ordered phase (3D XY). 



eluding the critical temperature and exponents. Com- 
bining these zero-temperature (T = 0) results with the 
existing results in the absence of quantum (temporal) 
fluctuations (K^ 1 = 0) or of the symmetry-breaking field 
(h = 0), we have constructed boundaries separating var- 
ious phases on the three planes (T = 0, K^ 1 = 0, and 
h = 0) in the three-dimensional (T, if ft,) space. The 
boundaries near T = and near K^ 1 = can then be 
speculated through the use of finite-size scaling and semi- 
classical methods, respectively, thus giving a schematic 
phase diagram in the (T, K -1 , h) space. 

We have also found via Monte Carlo simulations 
that the system with nonzero misfit undergoes a 
commensurate-incommensurate transition. The inter- 
chain interactions give rise to the correlations between 
neighboring CDWs in such a way that either the com- 
mensurate state or the incommensurate state is favored 
depending on the initial configuration: In the cooling pro- 
cess the CDWs remain incommensurate down to almost 



zero temperature while in the heating process a substan- 
tial amount of the commensurate CDWs survive thermal 
fluctuations. 

At strong thermal or quantum fluctuations, i.e., at high 
T or small K, the system is in the disordered phase. 
No CDW is formed and the system is expected to be 
metallic. On the contrary, weak fluctuations (low T and 
large K) favor the M-state clock ordered phase, in which 
commensurate CDWs are developed as long as 5 < 5 C . 
Accordingly, the interaction between the periodicity of 
the CDW and the underlying lattice periodicity drives 
the collective excitation to develop a gap, and the sys- 
tem becomes insulating. At moderate fluctuations or for 
large misfit (6 > S c ), on the other hand, incommensu- 
rate CDWs emerge. In this case the system may re- 
main conducting through collective Frohlich conduction, 
i.e., via sliding of the CDWs. Usually, the conductiv- 
ity via such collective modes is lower than that via un- 
condensated electrons in the disordered phase (without 
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FIG. 9: Schematic phase diagram for the coupled CDW sys- 
tem in the (T,t/-\<5) space, with the effects of the misfit 
taken into account. The surfaces depicted by thick solid lines 
and dashed lines separate the incommensurate CDW phase 
(IC) from the disordered phase (D) and from the commensu- 
rate CDW phase (C), respectively. 

CDW). In particular the CDW may be pinned in the 
presence of impurities, sharply decreasing the conductiv- 
ity. As the temperature is lowered, therefore, the sys- 



tem for weak quantum or thermal fluctuations becomes 
insulating via three possible routes, depending on the 
misfit and the inter-chain interaction: First, the com- 
mensurate CDW phase emerges directly from the high- 
temperature disordered phase; second, only the incom- 
mensurate CDW phase appears, reducing the conductiv- 
ity; third, the incommensurate CDW phase appearing 
first is followed by the commensurate phase emerging via 
the commensurate-incommensurate transition. 

Note that beginning with the Hamiltonian in Eq. 
we have taken into account only phase fluctuations and 
disregarded amplitude fluctuations. The latter are in 
general irrelevant in the RG sense, expected not to af- 
fect nature of the phase transition. On the other hand, 
there still lacks conclusive understanding of the dynamic 
properties in various phases. It is thus desirable to con- 
sider the responses to external electromagnetic fields, and 
for example, compute the conductivity, which can be ob- 
tained from the current or the average momentum in the 
presence of appropriate misfit. Detailed investigation of 
such dynamic responses is left for further study. 

Acknowledgments 

We thank G.S. Jeon for helpful discussions and ac- 
knowledge the partial support from the Korea Science 
and Engineering Foundation through the SKOREA Pro- 
gram and from the Ministry of Education of Korea 
through the BK21 Program. 



1 S. Kagoshima, H. Nagasawa, and T. Sambongi, One- 
Dimensional Conductors (Springer- Ver lag, Berlin, 1987). 

2 See, e.g., B. Grevin, Y. Berthier, G. Collin, and P. Mendels, 
Phys. Rev. Lett. 80, 2405 (1998); H.W. Yeom, et al., ibid. 
82, 4898 (1999); Y. Li, et al., ibid. 83, 3514 (1999). 

3 S. Coutinho, P. Pitanga, and P. Lederer, Phys. Rev. B 23, 
4567 (1981). 

4 G. Griiner, Density Waves in Solids (Addison- Wesley, 
Reading, 1994). 

5 M. Wallin, E.S. S0rensen, S.M. Girvin, and A.P. Young, 
Phys. Rev. B 49, 12115 (1994); S.L. Sondhi, S.M. Girvin, 
J.P. Carini, and D. Shahar, Rev. Mod. Phys. 69, 315 
(1997). 

6 S. Doniach, Phys. Rev. B 24, 5063 (1981); E. Simanek, 
ibid. 32, 500 (1985); J.V. Jose, ibid. 29, 2836 (1984); S. 
Kim and M.Y. Choi, ibid. 41, 111 (1990). See also A. van 
Otterlo, K.-H. Wagenblast, R. Fazio, and G. Schon, Phys. 
Rev. B 48, 3316 (1993); B.J. Kim, J. Kim, S.Y. Park, and 
M.Y. Choi, ibid. 56, 395 (1997). 

7 L. Onsager, Nuovo Cimento Suppl. 6, 279 (1949); R.P. 
Feynman, in Progress in Low Temperature Physics Vol. 1, 
edited by C. Gorter (North-Holland, Amsterdam, 1955), 
p. 17. 

8 G.A. Williams, Phys. Rev. Lett. 59, 1926 (1987); 61, 
1142(E) (1988). 

9 SR. Shenoy, Phys. Rev. B 40, 5056 (1989). 



10 G. Kohring, R.E. Shrock, and P. Wills, Phys. Rev. Lett. 
57, 1358 (1986); G. Kohring and R.E. Shrock, Nucl. Phys. 
B 288, 397 (1987); H. Kleinert, Gauge Fields in Condensed 
Matter Vol. I (World Scientific, Singapore, 1989), pp. 528- 
530; A.P. Gottlob and M. Hasenbusch, Physica A 201, 593 
(1993). 

11 J.M. Kosterlitz and D.J. Thouless, J. Phys. C 6, 1181 
(1973); J. Kosterlitz, ibid. 7, 1046 (1974). 

12 J.V. Jose, L.P. Kadanoff, S. Kirkpatrick, and DR. Nelson, 
Phys. Rev. B 16, 1217 (1977). 

13 K.G. Wilson and M.E. Fisher, Phys. Rev. Lett. 28, 240 
(1972). 

14 J. Cardy, Scaling and Renormalization in Statistical 
Physics (Cambridge University Press, Cambridge, 1996), 
pp. 100-101. 

15 M.E.J Newman and G.T. Barkema, Monte Carlo Methods 
in Statistical Physics (Oxford University Press, New York, 
1999), pp. 219-228. 

16 D.J. Scalapino, R.L. Sugar, and W.D. Toussaint, Phys. 
Rev. B 34, 6367 (1986). 

17 J. Rudnick, J. Phys. A 8, 1125 (1975); D. Kim and 
R.J. Joseph, ibid. 8, 891 (1975); T.W. Burkhardt, H.J.F. 
Knops, and M. denNijs, ibid. 9, L179 (1976); H.W.J. Blote 
and R.H. Swendsen, Phys. Rev. Lett. 43, 799 (1979); B. 
Nienhuis, E.K. Riedel, and M. Schick, Phys. Rev. B 23, 
6055 (1981); J. Lee and J.M. Kosterlitz, ibid. 43, 1268 



14 



(1991); P.D. Scholten and L.J. Irakliotis, ibid. 48, 1291 
(1993). 

W.Y. Shih and D. Stroud, Phys. Rev. B 28, 6575 (1983). 

See, e.g., M.Y. Choi, Physics of Complex Low-Dimensional 

Systems (Hanul, Seoul, 2000), pp. 49-57. 

P.S. Sahni, Phys. Rev. B 23, 1325 (1981). 

B. Topic, U. Haeberlen, and R. Blinc, Phys. Rev. B 42, 

7790 (1990). 



R.A. Craven and S.F. Meyer, Phys. Rev. B 16, 4583 
(1977); R.M. Fleming, D.E. Moncton, D.B. McWhan, and 
F.J. DiSalvo, Phys. Rev. Lett. 45, 576 (1981); J.-G. Yoon 
and S.-I. Kwun, Phys. Rev. B 35, 8591 (1987). 
Here we disregard the possible existence of re-entrance, 
which has been a source of controversy. See Ref. 0. 



